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INTRODUCTION 


To  most  accurately  approximate  any  physical  phenomenon,  a  mathematical  model  should 
be  posed  in  full  three-dimensional  space  and  time.  Even  with  the  advent  of  supercomputers, 
however,  obtaining  a  full  three-dimensional  solution  numerically  is  not  feasible  for  many 
investigators.  Access  to  a  supercomputing  facility  is  not  easily  obtainable  by  most.  When 
access  is  achieved,  time  and/or  memory  constraints  are  sometimes  imposed.  In  addition, 
software  advancements  are  still  lagging  behind  the  great  strides  being  made  in  hardware  design. 
This  makes  it  difficult  to  take  full  advantage  of  these  new  computer  architectures. 

Appropriately  reducing  the  dimensionality  of  the  mathematical  model  still  seems 
necessary.  The  best  scenario  would  be  a  reduction  to  one  spatial  dimension.  It  is  rare,  however, 
when  a  one-dimensional  problem  can  be  substituted  for  a  three-dimensional  one  without  an 
unacceptable  loss  in  physical  meaning.  As  a  result,  mathematical  models  in  two  spatial 
dimensions  are  still  predominant  in  most  scientific  fields. 

Solving  two-dimensional  problems  is  not  accomplished  without  some  level  of 
computational  difficulty.  Many  technological  situations,  modeled  in  two  spatial  dimensions  and 
time,  involve  the  rapid  formation,  propagation,  and/or  disintegration  of  small-scale  structures. 
Some  examples  are  shock  waves  in  compressible  flows,  shear  layers  in  laminar  and  turbulent 
flows,  phase  boundaries  in  nonequilibrium  processes,  combustion  fronts,  and  classical  boundary 
layers. 


Such  phenomena  often  arise  in  gun-related  problems  as  well.  For  example,  shock  waves 
occur  in  the  modeling  of  gas  dynamics  and  intemal/extemal  ballistics  while  swift  transients  are 
found  when  investigating  dynamic  effects  in  gun  tubes  (refs  1,2). 

Such  stmctures  pose  numerical  complications  since  standard  methods  require  locally  fine 
meshes  (in  space  and/or  time)  for  adequate  approximation  purposes.  Meanwhile,  the  location  of 
these  phenomena  is  generally  unknown  a  priori.  The  typical  response  is  a  globally  fine  mesh 
resulting  in  unnecessary  computational  effort  for  most  of  the  problem  domain. 

Adaptive  numerical  methods  were  developed  to  address  the  difficulties  associated  with 
problems  that  would  otherwise  require  dense  meshes  everywhere.  An  adaptive  method  is  one 
which  automatically  adjusts  its  solution  technique  so  that  an  answer  is  obtained  in  an  optimal 
manner.  Three  basic  adaptive  procedures  (/i-refinement,  r-refmement,  and  p-refmement)  have 
evolved.  All  have  proven  to  be  fairly  successful  at  optimizing  their  solution  processes,  although 
they  accomplish  this  task  through  different  means  (ref  3). 

To  the  uninitiated,  /i-refinement  (also  known  as  local  mesh  refinement)  and  adaptive 
methods  are  synonymous.  An  /i-refinement  approach  adjusts  the  spatial  and/or  temporal 
discretization  by  adding  and/or  deleting  mesh  points  so  that  finer  meshes  are  used  in 
neighborhoods  of  local  disturbances  while  coarser  meshes  are  used  elsewhere  (ref  4).  The  letter 
"h"  is  commonly  used  to  designate  the  distance  between  mesh  points,  hence  the  name. 
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An  r-refinement  procedure  adjusts  the  discretization  differently  than  an  /i-refmement 
scheme.  Instead  of  creating  new  points,  existing  points  are  moved  to  regions  of  high  activity 
from  regions  of  low  activity  as  time  progresses  (ref  5).  As  a  result,  r-refinement  is  also  known 
as  mesh  movement  for  time-dependent  problems.  In  static  problems,  mesh  points  are  said  to  be 
redistributed,  as  opposed  to  moved,  consequently  the  letter  "r"  is  used  to  describe  these  types  of 
algorithms. 

The  basic  philosophy  behind  p-refinement  is  different  fi'om  the  other  two  methods. 
Algorithms  based  on  /?-refinement  or  r-refinement  attempt  to  modify  the  level  of  discretization 
in  some  way.  The  identical  numerical  solver  is  used  throughout  the  entire  spatial  domain  and  for 
all  time. 

In  p-refinement  schemes,  the  idea  is  to  utilize  numerical  methods  with  varying  orders  of 
accuracy  in  different  subdomains,  as  opposed  to  directly  adjusting  mesh  points  and/or  time  steps. 
For  example,  a  finite  element  method  with  p-refinement  capabilities  would  automatically  decide 
what  degree  polynomial  (linear  (p  =  1),  quadratic  (p  =  2),  cubic  (p  =  3),  etc.)  would  be 
appropriate  to  use  as  basis  functions  in  various  regions  of  the  problem.  In  general,  low  order 
methods  would  be  used  in  regions  where  the  solution  is  changing  very  little  (i.e.,  where  solution 
gradients  are  small).  The  higher  order  methods  would  be  employed  when  and  where  solution 
characteristics  change  dramatically,  as  in  the  cases  described  above  (ref  6). 

All  three  refinement  schemes  have  their  individual  strengths  and  weaknesses  (ref  3).  For 
example,  for  a  fixed  level  of  discretization  and  a  fixed  order  of  accuracy,  an  r-refinement 
procedure  will  perform  optimally.  However,  once  this  optimization  is  achieved,  a  finer 
discretization  and/or  a  higher  order  method  must  be  employed  in  order  to  increase  accuracy.  As 
a  result,  these  adaptive  strategies  are  often  combined. 

An  adaptive  finite  element  method  that  performed  both  /z-refinement  and  r-refinement 
was  developed  for  time-dependent  problems  in  one  spatial  dimension  (ref  7).  This  code  has 
made  contributions  in  various  areas  related  to  cannon  development  such  as  gun  dynamics, 
impact  and  penetration,  and  rail  gun  technology  (ref  8).  However,  the  dimensionality  restriction 
has  proven  to  be  too  severe  for  this  method  to  evolve  into  a  practical  tool  for  gun  problems. 

The  goal  has  become  to  extend  this  work  to  two  spatial  dimensions.  In  particular,  a  two- 
dimensional  mesh  moving  technique  for  rectangular  domains  was  suggested  by  the  one¬ 
dimensional  results  (ref  9).  The  purpose  of  this  repon  is  to  examine  the  properties  and 
capabilities  of  this  two-dimensional  moving  scheme.  An  extension  to  other  geometries  is 
implied  by  the  two-dimensional  results,  but  for  simplicity,  this  report  will  deal  with  only 
rectangular  regions.  Reports  concerned  with  more  general  geometries  will  be  forthcoming. 

Ultimately,  other  adaptive  techniques  will  be  paired  with  r-refinement  in  order  to 
produce  a  more  robust  computational  tool.  Since  /i-refinement  and  p-refinement  normally 
involve  solving  a  given  problem  more  then  once  at  selected  times,  they  are  more  costly  than 
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r-refinement.  Proper  mesh  movement  can  postpone  the  necessity  for  additional  adaptivity  to 
take  place  and  hence  reduce  the  computational  costs. 

Significant  accomplishments  were  achieved  in  the  development  of  the  aforementioned 
one-dimensional  r-refinement  procedure  (ref  9).  First  and  foremost  was  the  mesh  moving 
algorithm  itself.  This  scheme  has  proven  to  be  robust,  reliable,  stable,  easily  controllable,  and  a 
natural  partner  with  other  refinement  techniques.  In  addition,  a  stability  criterion  was  derived 
demonstrating  why  some  other  schemes  become  unstable  and  illustrating  how  to  construct 
different  types  of  schemes  with  various  stability  properties. 

The  two-dimensional  r-refinement  procedure  analyzed  here  is  a  relatively  simple 
extension  of  the  one-dimensional  algorithm.  Due  to  the  nature  of  this  extension,  all  of  the  "nice" 
properties  of  the  one-dimensional  movement  (reliability,  stability,  controllability,  etc.)  are 
inherited  by  the  two-dimensional  scheme. 

In  the  next  section,  theoretical  aspects  of  the  two-dimensional  mesh  moving  algorithm 
are  presented.  First,  the  one-dimension^  algorithm  is  reviewed .  Then,  the  extension  to  two- 
dimensional  rectangular  domains  is  detailed.  Concepts  such  as  existence,  uniqueness,  initial 
condition  dependency,  and  stability  are  discussed  and  the  manner  in  which  the  two-dimensional 
scheme  inherits  these  properties  from  the  one-dimensional  scheme  is  outlined.  An  example  is 
provided  which  illustrates  the  stability  properties  of  this  method  as  compared  to  another  more 
obvious  extension  to  two  dimensions. 

The  third  section  entitled,  “Control,”  discusses  a  more  practical  aspect  of  the  two- 
dimensional  mesh  moving  scheme.  Ultimately,  the  scheme  must  be  coupled  with  a  numerical 
partial  differential  equation  (PDF)  solver  (ref  10).  In  one  dimension  it  was  discovered  that  mesh 
movement  could  degrade  the  PDE  solver’s  performance  by  moving  points  too  dramatically. 
Space-time  elements  could  become  too  distorted  for  the  numerical  solver  to  resolve  adequately. 
Some  control  over  mesh  movement  had  to  be  exercised  in  order  to  avoid  these  situations  (ref  7). 
The  same  problem  is  anticipated  in  two  dimensions,  and  this  section  describes  how  to  exercise 
control  over  the  two-dimensional  mesh  moving  algorithm  in  a  relatively  simple  manner.  First, 
space-time  distortion  is  addressed,  followed  by  a  discussion  on  the  resolution  of  space-space 
distortion.  Examples  illustrating  the  utility  of  this  control  methodology  are  presented. 

The  fourth  section,  “Discussion  and  Conclusions,”  reviews  the  results  of  the  previous 
sections  and  draws  conclusions.  Weaknesses  as  well  as  strengths  of  the  two-dimensional  mesh 
moving  scheme  are  discussed  and  future  plans  detailed.  This  section  is  followed  by  a  list  of 
references  and  a  sequence  of  figures,  which  are  referred  to  within  the  text  of  this  report. 

THEORY 

The  one-dimensional  mesh  moving  algorithm  is  based  on  the  notion  of  equidismbution 
(ref  11).  In  this  context,  a  one-dimensional  equidistribution  problem  is  the  determination  of  a 
dynamic  partition 
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n/r)  =  ia  =  Xq  <  X^it)  <  X^it)  <  ...  <  ;C/.i(r)  <  x,  =  b) 


(1) 


of  (a,b)  into  I  elements  such  that 


w^(x,t)dx  =  K(0  =  4  f  ,1  =  1, 2, ...,/,  r  ^  0 . 


(2a) 


or  equivalently 


r‘^‘^w^(x,t)dx  =  i  K(0  =  4  f  ^yvXx,t)dx ,  i  =  1, 2, ...,  /,  r  ^  0 , 

Ja  IJa 


(2b) 


where  the  weight  function  w,(x,t)  >0,x  e  [a, b],  t  ^0,  is  usually  dependent  on  the  solution  of 
the  underlying  PDE.  For  example,  Wj  has  been  chosen  to  be  proportional  to  the  solution's 
gradient,  curvature,  and  local  spatial  discretization  error  (refs  12-20). 

Reducing  the  global  discretization  error  is  the  ultimate  goal,  and  towards  that  end 
equidistributing  the  local  discretization  error  would  be  the  best  procedure  to  follow.  However, 
estimates  of  the  total  discretization  error  are  not  easy  to  make,  even  locally,  and  this  has  led  to 
these  other  choices  for  the  weight  function.  Since  for  most  numerical  methods  the  discretization 
error  is  proportional  to  some  order  of  derivative  of  the  solution,  this  has  become  a  popular 
choice  (refs  13,14,15,16,17,20).  It  has  also  been  argued  that  for  certain  problems,  an  appropriate 
physical  characteristic  such  as  density  is  a  suitable  weight  function  (ref  18).  When  a  method  of 
lines  approach  is  employed,  the  local  spatial  discretization  error  can  be  used  since  the  ordinary 
differential  equation  solver  can  reduce  the  temporal  error  component  to  a  relatively  insignificant 
level  (refs  12,19). 


It  was  discovered  that  Eq.  (2)  was  subject  to  unstable  behavior  as  written  (ref  21).  In 
order  to  construct  a  stable  mesh  moving  scheme,  it  was  necessary  to  define 


f  Wj(^,r)flf^ 
J  a 


(3) 


rewrite  Eq.  (2b)  as 


<^(x.(t),t)  -1=0,1  =  0,1,2,...,/, 


(4) 


and  differentiate  Eq.  (4)  with  respect  to  time  to  obtain 
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=  0  ,i  =  0,1,2,...,/, 
dt 


(5a) 


a>(x.(0),o)  =  0° ,  /  =  0,1,2, 


(5b) 


With  emphasis  placed  on  the  stability  of  the  functions  ^(x/t),t),  I  =  it  is  seen  that 

Eq.  (5)  is  neutr^ly  stable  (ref  9).  Initial  perturbations  of  6,  (through  a:,),  I  =  will  neither 

grow  nor  decay. 

Another  advantage  to  the  formulation  of  Eq.  (5),  is  that  as  a  result  of  the  differentiation, 
an  initially  equidistributed  mesh  is  no  longer  necessary.  Any  initial  partition  determines 
Eq.  (5b)  uniquely,  and  Eq.  (5a)  guarantees  that  the  resulting  movement  will  be  stable  and 
optimal,  although  the  resulting  mesh  will  not  necessarily  be  equidistributed.  Equidistribution  is 
guaranteed  for  all  time,  however,  if  the  initial  mesh  is  an  equidistributed  one. 

As  mentioned  previously,  extending  the  neutrally  stable  mesh  moving  scheme  (cf., 

Eq.  (5))  to  two-dimensional  space  can  be  accomplished  with  all  of  the  stability  properties  left 
intact  (ref  9).  Consider  a  rectangular  two-dimensional  domain 

Q  =  {(a,}’)  \  a<x<b,c^y^d}  (6) 


and  the  extension  of  Eq.  (5)  to  Q.  To  that  end,  define 

f''f^W2(^,r],t)dr]d^ 

Wixyt)  =  IS-lS - 

f‘’f\a,r]d)dr]d^ 

J  fl  J  C 


(7) 


where  W2{x,y,t)  is  a  positive  weight  function  on  Q  such  that  T(A,y,rj  is  a  well-defined  function 
with  the  necessary  continuity  properties.  Furthermore,  let 

=  { (A.(0,y/0)  i  /  =  0, 1, ...,/  ,  ;■  =  0, 1, ...,/}  (8) 

denote  a  partition  of  Q  into  IxJ  subrectangles  at  any  time  t  such  that 
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a  =  Xq<  Xj(0  <  ...  <  x,_^{t)  <Xj  =  b  , 

(9a) 

c  =  yo<  y:(^)  <  -  <  <yj  =  d . 

(9b) 

In  order  for  n;y(0  to  move  in  a  neutrally  stable  manner  for  any  initial  partition  Xly/O), 
demand  that  it  obey  the  following  equations  for  r  >  0: 

4mx,(r),d,r)]  =  0  ,  /  =  1,2 . 7-1  , 

dt 

(10a) 

■^mb,yp),t)]  =  0  ,y  =  1,2,. ..,7-1  , 
dt  ^ 

(10b) 

^(x.(0),d,0)  =  X.  ,  i  =  1,2,...,/-!  , 

(10c) 

m^/Ol.O)  =  Y.J  =  1,2,..., 7-1  . 

(lOd) 

The  existence,  uniqueness,  and  stability  of  the  two-dimensional  movement  is  guaranteed  by  the 
one-dimensional  theory.  By  substituting 

w^(x,t)  =  j'‘w2(x,r],t)dr\ 

(11a) 

into  Eq.  (5),  it  can  be  seen  that  Eq.  (10a)  is  equivalent  to  neutrally  stable,  one-dimensional 
movement  in  the  x-direction.  Similarly,  by  replacing  the  variable  x  with  the  variable  y  and  the 
interval  (a,b)  with  the  interval  {c,d)  in  Eq.  (5),  and  substituting 

H'i(y,0  =  (lib) 

a 


it  can  be  seen  that  Eq.  (10b)  is  equivalent  to  neutrally  stable,  one-dimensional  movement  in  the 
y-direction. 

As  in  the  one-dimensional  case  (ref  9),  the  normalization  procedure  (cf.,  Eq.  (7))  is 
crucial  for  stability  in  the  construction  of  the  two-dimensional  moving  scheme,  Eq.  (10).  If  such 
a  normalization  is  not  performed  the  resulting  equations  are 
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4[0(^,W.rf,O]  =  ,  i  =  1,2,...,/-!  ,  (12a) 

dt  I  dt  ^ 

4-\.^{b,yft),t)]  =  {4-^Q{b4,t)]  ,j  =  1,2,...,  7-1  ,  (12b) 

dt  ^  J  dt 

0(x.(O),d.O)  =  j0(b,d.O)  ,  i  =  1,2 . 7-1  ,  (12c) 

&(b,yj(0),0)  =  ^©(M,0)  ,y  =  1,2,...,  7-1  ,  (12d) 


where 


f^W2(i,ri,t)drid^  . 

J  a  c 


(12e) 


Such  a  system  is  not  stable  when  W2(x,y,t)  is  a  decreasing  function  of  time  as  illustrated  by  the 
following  example. 

Examp.l£..l 

Consider  the  two-dimensional  heat  equation 

0<a:<1,  0<y<l,  t>0 

subject  to  the  initial  condition 

u(x,y,0)  =  sinTtjc  sinTty,  0<x<l,  0<y<l 

and  boundary  conditions 

u(0,y,t)  =  u(l,y,t)  =  0,  0<y<l, 
u(x,0,t)  =  u(Xyl,t)  =  0,  0<a:<1,  tiO 


(13a) 

(13b) 


(13c,d,e,f) 


The  exact  solution  of  this  problem  is 
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(14) 


71^ 

u(x,y,t)  =  sm(7ij:)  sinCu}^)  exp(- — r),  0<a:<1,  0<j<1 

4 

For  the  purposes  of  mesh  movement  take 

W2ix,y,t)  =  uix,y,t).  (15) 

Since  this  problem  and  W2(x,y,t)  are  separable,  the  correct  strategy  is  to  generate  an 
equidistributed  mesh  at  time  t  =  0  and  use  it  for  all  time.  Hotvever,  W2(x,y,t)  is  a  decreasing 
function  of  time  and,  thus,  the  solutions  of  Eq.  (12)  are  expected  to  be  unstable  (ref  21). 

For  the  given  Wi,  exact  solutions  to  Eqs.  (10)  and  (12)  can  be  found.  The  exact  solution 
for  Eq.  (10)  is 

x.(r)  =  a:.(0),  yff)  =  y^.(0),  i  =  1, 2,  ...,/-l,  y  =  1, 2,  ...,/-l  (16a, b) 

The  exact  solution  for  Eq.  (12a)  is 

x-{t)  =  ^arccos[a.-  (a.-  cosTtxf)  exp(-jr)],  i  =  1, 2,  ...,/-l  (17a) 

where 

a.  =  1  -  2-y  and  xf  =  x.(0),  i  =  1, 2,  ...,/-l .  (17b, c) 


The  solution  to  Eq.  (12b)  can  be  found  by  substituting  y  for  x,  j  for  I,  and  J  for  7  in  Eq.  (17). 

The  trajectories  for  the  meshes  produced  by  Eqs.  (16)  and  (17)  are  displayed  in  Figures  1  and  2, 
respectively,  and  the  unstable  behavior  of  Eq.  (17)  is  clearly  visible.  A  mesh  with  values  of 
1  =  J  =  3  was  initially  equidistributed,  and  initial  perturbations  of  0.0]  for  7  =  y  =  7  and  of  -O.OI 
for  7  =  y  =  2  were  introduced. 

It  should  be  noted  that  Eq.  (10)  could  easily  be  adapted  to  nonrectangular  regions  if 
another  natural  set  of  coordinates  existed  for  the  given  problem  domain,  e.g.,  polar  coordinates 
for  circular  regions  or  boundary-fitted  coordinates  as  used  for  flow  around  air  foils  (ref  22). 

This  natural  coordinate  system  could  be  used  for  the  mesh  equations  (cf.,  Eq.  (10))  as  well.  All 
that  is  required  is  to  transform  Eq.  (10)  from  rectangular  coordinates  to  the  more  natural  system. 
(This  is  not  the  extension  to  nonrectangular  geometries  alluded  to  in  the  “Introduction,”  but 
rather  an  extension  to  nonrectangular  coordinate  systems  for  appropriately  defined  regions.) 
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CONTROL 


In  the  one-dimensional  situation,  a  positive  constant  can  be  added  to  Wy  to  ensure  the 
existence  of  a  unique  solution.  It  should  be  noted  that  adding  a  positive  constant  to  will 
ensure  the  existence  of  a  unique  solution  to  Eq.  (10)  as  well.  It  was  discovered  that  the  inclusion 
of  this  constant  also  provided  a  means  in  which  to  control  the  one-dimensional  movement 
(ref  8). 


Control  was  deemed  necessary  since  mesh  movement,  as  defined  by  Eq.  (5),  could  be  too 
dramatic  and  amplify  the  discretization  error  (especially  the  temporal  component)  of  a  numerical 
PDE  solver  by  overly  distorting  the  mesh  (ref  7).  In  the  one-dimensional  case,  there  was  only 
space-time  distortion  to  be  concerned  about.  Now  in  two  dimensions  there  exists  the  possibility 
of  space-space  distortion  as  well.  (Since  in  both  one-  and  two-dimensional  space  there  is  the 
well-recognized  problem  of  properly  transitioning  from  a  dense  mesh  to  a  sparse  one,  it  will  not 
be  discussed  here.)  External  control  is  even  more  of  a  necessity  in  two  dimensions,  and  the 
addition  of  a  positive  constant  to  W2  appears  to  be  a  simple  method  of  accomplishing  the  task. 

For  one-dimensional  mesh  movement,  it  was  discovered  that  the  solution  to  Eq.  (5),  with 
a>  0  added  to  Wj,  could  be  written  as  a  linear  combination  of  two  other  solutions  to  Eq.  (5)  with 
appropriate  weight  functions  (ref  8).  This  is  also  the  case  in  two  dimensions  since  Eq.  (10),  as 
previously  stated,  represents  one-dimensional  movement  in  the  respective  coordinate  directions. 
In  order  to  demonstrate  this  fact,  consider  the  jc-component  of  the  two-dimensional  moving 
scheme  (cf.,  Eq.  (10a)).  The  y-component  behaves  similarly. 


Letx/r),  qi(t),  and  Ui,  I  =  0, 1 , I,  denote  mesh  trajectories  for  1 2:  0  such  that 

i  r  ^  r 

'a  J  c 


M[w2(^,Ti,r)  +  a]dy\dl  =  -  +  a]drid^ ,  /  =  0, 1, ...,/ ,  (I9a) 

J  a  J  c  I  J  a  J  c 


f‘'<^U^W2a,ri,t)dr\dl  =  i  =  0,l, (I9b) 

J  a  c  I  J  a  J  c 


I  r  b  r  d 


and 


f‘^adr]d^  =  1  T*  f^adr]d^  ,  i  =  0, 1, 

J  a  J  c  I J  o  c 


(19c) 


Note  that  the  trajectories,  I  =  0,1, ...,  /,  are  solutions  to  the  equidistribution  problem  as 
originally  defined  and  w.,  /  =  0,  i, ...,/,  are  constant  trajectories  defining  a  uniform  mesh. 

The  solutions,  x/t),  7  =  0,  7, ...,/,  are  actually  linear  combinations  of  the  equidistributed 
trajectories,  q^t),  1  =  0,1, ....  7,  and  the  constant  trajectories,  I  =0,1, ...,  7 ,  according  to  the 
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one-dimensional  theory.  Furthermore,  the  dominant  component  is  determined  by  the  following 
relationship  (ref  8): 


\x^it) 


max  f^W2(x,T],t)dr\ 


U.\ 


a  (d  -  c) 


qft)  ,  /  =  0, 1,...,/. 


(20) 


For  large  a,  as  compared  to  the  maximum  value  of  the  solutions,  x/r),  I  =  0,1, ...,  /,  are  very 
close  to  being  constant  trajectories.  Conversely,  for  small  a,  the  solutions,  I  =  0,1, ....  /, 
are  very  close  to  the  equidistributed  positions.  (N.B.,  for  a  =  0,  Eq.  (19a)  and  Eq.  (19b)  are 
identical.)  Appropriate  choices  of  a  determine  the  amount  of  movement  deemed  necessary  as 
illustrated  by  the  following  example. 

Example  2 

Consider  the  first-order  wave  equation 

M,  +  piu^  +  Uy)  =  0,  0<x<l,  0<y<l,  t>0,  (21a) 

with  initial  and  boundary  conditions  and  the  constant,  p,  selected  so  that  the  exact  solution  is 

u(x,y,t)  =  ^  [1  “  tanh(10x  +  lOy  -  20t)],  (2lb) 


which  is  an  oblique  wave  front  (slope  =  -1)  that  is  initially  centered  on  the  line 

lOx  +  lOy  =  0  (21c) 


and  subsequently  moves  across  the  domain,  =  [0,l]x[0,l],  from  left  to  right,  i.e.,  in  the 
direction  of  increasing  x.  If  one  wished  to  solve  this  problem  numerically,  incorporating  mesh 
movement,  a  reasonable  choice  for  the  weight  function,  W2,  would  be: 

H'2(x,y,0  =  sech^(10x  +  lOy  -  200  +  a  «  uf;x,y,t)  +  a.  (22) 

Eq.  (10)  was  solved  using  the  differential-algebraic  solver,  DASSL  (ref  23),  with  the 
above  choice,  Eq.  (22),  for  the  weight  function,  Wj,  and  three  different  values  for  the  constant,  a. 
In  all  instances,  the  initial  condition  was  a  uniform  partition  with  7  =  7=7.  The  purpose  was  to 
illustrate  the  effect  the  parameter  a  has  on  mesh  moving  dynamics,  as  well  as  the  amount  of 
control  available  to  the  user. 

Meshes  at  time,  t  =  0.25,  for  values  of  a  =  0.02, 0.1,  and  0.5  are  displayed  in  Figures  3, 
4,  and  5,  respectively.  The  gray  area  of  the  figures  represents  the  region  of  the  domain  behind 
the  front,  while  the  region  ahead  of  the  front  remains  white.  As  predicted,  movement  from  the 
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initially  uniform  positions  is  more  dramatic  for  a  =  0.02,  while  for  a  =  0.5  mesh  points  have 
barely  moved  at  all. 

The  intermediate  value  of  a  =  0.1  would  seem  to  be  optimal.  Mesh  points  have  moved 
significantly  so  as  to  address  the  issue  of  the  moving  wave  front  adequately.  On  the  other  hand, 
the  movement  is  not  so  dramatic  so  as  to  possibly  disrupt  the  ability  of  a  PDE  solver  by  causing 
excessive  space-time  distortion.  In  general,  one  would  expect  intermediate  values  of  a  to  be  best 
with  specific  problems  and  solvers  demanding  more  extreme  values  on  occasion. 

In  the  above  example,  spatial  distortion  of  the  computational  cells  was  not  discussed. 
Since  the  moving  front  intersected  both  the  x-  and  y-axes  at  45  degrees,  the  aspect  ratio  of  the 
cells  for  a  =  0.1  and  0.5  did  not  change  appreciatively.  For  a  =  0.02,  mesh  movement  tended  to 
elongate  the  cells  in  one  direction  or  the  other.  It  is  conceivable  that  this  space-space  distortion 
could  become  great  enough  so  as  to  affect  the  accuracy  of  a  PDE  solver. 

In  order  to  address  this  problem,  it  should  be  noted  that  there  is  no  requirement  that  a  be 
the  same  value  in  Eq.  (10a)  as  it  is  in  Eq.  (10b).  (For  that  matter,  there  is  no  requirement  that  the 
entire  weight  function  be  the  same  in  the  two  equations.)  Different  values  for  the  constant,  a, 
can  be  used  to  deal  with  spatial  mesh  distortion,  if  necessary.  This  fact  is  illustrated  by  the 
following  example. 

Example  3 

Once  again,  consider  the  first-order  wave  equation 

+  p(m^  +  up  =  0,  0<j:<1,  0<y<l,  r>0,  (23a) 

but  with  initial  and  boundary  conditions  and  the  constant,  p,  selected  so  that  the  exact  solution  is 

u{x,y,t)  -  ~  tanh(10j:  +  2y  -  lOr)],  (23b) 

which  is  an  oblique  wave  front  (slope  =  -5)  that  is  initially  centered  on  the  line 

lOx  +  2y  =  0  (23c) 


and  subsequently  moves  across  the  domain,  Q  =  [0,l]x[0,l],  from  left  to  right,  i.e.,  in  the 
direction  of  increasing  x.  If  one  wished  to  solve  this  problem  numerically,  incorporating  mesh 
movement,  a  reasonable  choice  for  the  weight  function,  W2,  would  be; 

W2ix,y,t)  =  sech^(10x  +  2y  -  lOr)  +  a  «  ufx^y,t)  +  a.  (24) 
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Eq.  (10)  was  solved  using  the  differential- algebraic  solver,  DASSL  (ref  23),  with  the 
above  choice,  Eq.  (24),  for  the  weight  function,  with  a  =  0.1  (in  both  equations),  and  for  an 
initially  uniform  mesh  with  /  =  /  =  7.  The  resulting  mesh  at  time,  r  =  0.5,  is  displayed  in 
Figure  6.  The  gray  area  of  the  figure  represents  the  region  of  the  domain  behind  the  front,  while 
the  region  ahead  of  the  front  remains  white.  The  elongation  of  the  computational  cells  is  clearly 
evident. 


Eq.  (10)  was  then  solved  using  Eq.  (24)  and  an  initially  uniform  mesh  with  I  =  J  =  7, 
but  with  a  =  1.0  in  Eq.  (10a)  and  a  =  0.1  in  Eq.  (10b).  The  resulting  mesh  at  time,  t  =  OJ,  is 
displayed  in  Figure  7.  The  gray  area  in  the  figure  represents  the  region  of  the  domain  behind  the 
front,  while  the  region  ahead  of  the  front  remains  white.  The  aspect  ratio  of  the  computational 
cells  has  not  been  altered  significantly  from  the  initial  aspect  ratio  of  1. 


DISCUSSION  AND  CONCLUSIONS 

A  neutrally  stable  two-dimensional  mesh  moving  procedure  was  proposed  (cf,  Eq.  (10)). 
Based  on  a  successful  one-dimensional  algorithm  (cf ,  Eq.(5)),  this  procedure  can  be  easily 
implemented  in  a  two-dimensional  problem  domain  for  which  a  natural  set  of  coordinates  exist, 
e.g.,  rectangular  domains  and  Cartesian  coordinates.  Theoretical  issues  (existence,  uniqueness, 
etc.)  were  discussed  and  it  was  demonstrated  how  these  properties  are  inherited  from  the  one¬ 
dimensional  formulation.  Examples  were  presented  illustrating  the  two-dimensional  mesh 
movement’s  attributes  of  stability  and  controllability.  A  more  general  extension  to  other 
geometries,  regardless  of  the  coordinate  system,  has  been  formulated  but  this  report  concentrated 
on  the  extension  to  rectangular  regions  in  order  to  introduce  the  basic  concepts. 

A  mesh  moving  scheme  based  on  Eq.  (10)  would  have  various  strengths  and  weaknesses 
just  like  any  numerical  method,  adaptive  or  otherwise.  For  example,  such  a  scheme  would  be 
extremely  stable  and  could  be  put  under  user  control  with  minimi  interface  requirements.  Any 
initial  mesh  would  be  a  valid  initial  condition  and  mesh  points  would  be  guaranteed  not  to 
coalesce.  Unacceptable  mesh  distortion  could  occur,  but  both  space-time  and  space-space 
distortion  could  easily  be  addressed  by  the  user. 

Individual  point  movement  would  not  be  independent  in  the  strictest  sense  of  the  word, 
however.  A  routine  examination  of  Eq.  (10)  reveals  that  Eq.  (10a)  is  not  dependent  on  y,  while 
Eq.  (10b)  has  no  x  dependency.  This  implies  that  all  points  on  the  same  coordinate  line  will 
move  in  unison  in  the  opposite  coordinate  direction  (e.g.,  all  points  with  identical  x-component 
values  will  generate  velocity  and  displacement  vectors  with  identical  y-component  values,  and 
vice  versa).  In  a  sense  then,  Eq.  (10)  delineates  a  line-moving  scheme  rather  than  a  point- 
moving  one.  This  makes  it  difficult  for  such  a  method  to  align  points  properly  when 
encountering  structures  oblique  to  the  coordinate  axes.  Orienting  the  initial  mesh  so  as  to  align 
properly  with  the  geometry  of  some  known  structure  can  overcome  this  weakness,  but  it  is  rare 
that  such  detailed  information  is  initially  available. 
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A  mesh  moving  scheme  with  individual  and  independent  point  movement  has  been 
formulated.  This  independent  scheme  is  also  neutrally  stable,  but  more  susceptible  to  mesh 
distortion.  A  thorough  analysis  of  this  independent  scheme  and  a  comparison  with  the  algorithm 
described  in  this  report  is  forthcoming.  Included  with  this  analysis  will  be  the  development  of  a 
two-dimensional  stability  criterion  analogous  to  the  one-dimensional  stability  criterion  for  mesh 
movement  based  on  equidistribution.  The  stability  analysis  presented  here  is  specific  to  the 
particular  two-dimensional  moving  scheme  described  herein. 

Issues  such  as  how  best  to  evaluate  the  weight  function,  and  how  to  properly  couple 
this  moving  scheme  with  a  PDE  solver  have  not  been  addressed  in  this  report.  Answers  to  such 
questions  are  dependent  on  the  properties  and  characteristics  of  the  PDE  solver  involved.  It  is 
best  to  delay  dealing  with  those  issues  until  after  a  particular  PDE  solver  has  been  chosen.  In 
the  future,  it  is  planned  to  combine  this  mesh  moving  scheme  with  some  existing  PDE  solver 
and  these  issues  will  be  reported  on  then,  as  well  as  the  relative  success  or  failure  of  such  a 
coupling. 
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Figure  1  Stable  mesh  trajectories  for  Example  1 
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Figure  2  Unstable  mesh  trajectories  for  Example  1 


Figure  3  Mesh  for  Example  2  at  f  =  0.25  and  a  =  0.02 
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Figure  6  Mesh  for  Example  3  at  f  =  0.5  and  a  *  0.1  in  both  Eqs.  (10a, b) 
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